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Abstract 


The study combines available data on species distribution in eastern Indochina to investigate the phylogeographical genetic and mor- 
phological diversity of tree rats (Chiromyscus, Rodentia, Muridae) and to specify their natural ranges. We examined the diversity and 
distribution of tree rats over its range, based on recent molecular data for mitochondrial (Cyt b, CO/) and nuclear (JRBP, RAG/ and 
GHR) genes. The study presents the most complete and up-to-date data on the distribution and phylogeography of the genus in east- 
ern Indochina. As revealed by mitochondrial genes, C. /angbianis splits into at least four coherent geographically-distributed clades, 
whereas C. thomasi and C. chiropus form two distinctive mitochondrial clades each. Chiromyscus langbianis and C. chiropus show 
significant inconsistency in nuclear genes, whereas C. thomasi shows the same segregation pattern as can be traced by mitochondrial 
markers. The Northern and Southern phylogroups of C. thomasi appear to be distributed sympatrically with northern phylogroups of 
C. langbianis in most parts of eastern Indochina. The mitochondrial clades discovered are geographically subdivided and divergent 
enough to suspect independent subspecies within C. langbianis and C. thomasi. However, due to the insufficiency of obvious mor- 
phological traits, a formal description is not carried out here. The processes of recent fauna formation, species distribution patterns, 
dispersion routes and possible natural history in Indochina are discussed. 
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Introduction (Thomas, 1891), C. langbianis and C. thomasi Balakirev, 


Abramov & Rozhnov, 2014. 


Genus Chiromyscus Thomas, 1925 is currently assigned 
to the Dacnomys division of the tribe Rattini (Musser and 
Carleton 1993, 2005). It was first described, based on Mus 
chiropus (Thomas, 1891) from East Burma (= Myanmar) 
and for a long time was considered monotypic. Its close 
relationships with Niviventer langbianis (Robinson, 
Kloss, 1922) and N. cremoriventer (Miller, 1900) were 
initially suspected and discussed by Musser (Musser 
1973, 1981; Musser and Carleton 1993, 2005). The genus 
was recently re-established by Balakirev et al. (2014) 
as comprising at least three recent species: C. chiropus 


To date, Thomas’ tree rat, C. thomasi, 1s known to be 
distributed from southwest China (Lan et al. 1994; Chen 
et al. 1995; Wang 2002) to eastern Myanmar and north- 
ern Thailand (Marshall 1977), Vietnam (Dang et al. 1994; 
Lunde and Nguyen 2001; Can et al. 2008, Balakirev et 
al. 2014) and Laos (Musser 1981; Corbet and Hill 1992; 
Aplin et al. 2003; Aplin and Lunde 2016), where it has 
long been known under the name C. chiropus. The Da- 
lat tree rat C. langbianis was described from the Dalat 
Plateau in southern Vietnam (Robinson and Kloss 1922) 
and is currently recorded throughout Vietnam and Laos 
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(Musser 1973; Dang et al. 1994; Lunde et al. 2003; Bal- 
akirev et al. 2011; Lu et al. 2015) and southern China 
(Lu et al. 2015), including Hainan Island (Ge et al. 2018). 
Fea’s tree rat C. chiropus 1s known from Myanmar and 
easternmost India (Musser 1973) and southern Vietnam 
(Balakirev et al. 2014). 

The evolution of the genus Chiromyscus was affected 
by the natural history of this region. Continuous forest 
cover in Indochina existed during a considerable pro- 
portion of the Pliocene and Pleistocene (Metjaard and 
Groves 2006), enabling the direct contact, dispersion 
and genetic exchange of western and eastern Indochinese 
populations, which are currently mostly interrupted by 
the extensive deforestation in areas of central Thailand 
and southern Cambodia. During Pleistocene glaciations, 
the forest edge was located at lower elevations (Cox and 
Moore 2000) and substantial forest contraction happened 
during the last glacial periods, as evidenced in Peninsular 
Malaysia and Palawan (Wurster et al. 2010). Cranbrook 
(2000) and Gathorne-Hardy et al. (2002) stressed that 
most areas in this region were covered with savannah 
vegetation unfavourable to arboreal species during long 
periods of the Quaternary epoch. However, the most re- 
cent surveys on the biogeography and paleoenvironments 
of the Sunda Shelf have suggested that the interactions 
between climate and sea level and their effects on the dis- 
tribution patterns of the fauna and flora are more complex 
(van den Bergh et al. 2001; Metjaard and Groves 2006) 
and environmental conditions changed repeatedly during 
the Pleistocene. These global processes of ecosystems 
change were likely to have had an impact on the recent 
genetic structure of Chiromyscus species. This group of 
rodents 1s still quite rare in museum collections, so little 
information is available about its natural diversity, dis- 
tribution and phylogenetic relationships and only a few 
specimens have been genetically characterised. In this pa- 
per, we combine available data, including novel data, on 
species distribution in eastern Indochina to investigate the 
phylogeography and diversity of these species, both ge- 
netic and morphological and specify their natural ranges. 


Material and methods 


Specimens and samples 


A great number of Chiromyscus specimens, obtained in 
Vietnam during a series of field expeditions of the Joint 
Russian-Vietnamese Tropical Research and Technologi- 
cal Centre between 2007 and 2018, were sampled for ge- 
netic analysis in full agreement with current Vietnamese 
regulations in the field of Nature Protection and Biodi- 
versity Conservation. We followed the guidelines of the 
American Society of Mammalogists during the collection 
and handling of the animals used in this survey (Gannon 
et al. 2011). The museum specimens were kept in the Zo- 
ological Museum, Moscow State University, Moscow, 
Russia (ZMMU) and the Zoological Institute, Russian 
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Academy of Sciences, Saint-Petersburg, Russia (ZIN); 
genetic samples are part of collection of Joint Russian— 
Vietnamese Tropical Research and Technological Centre, 
Hanoi, Vietnam. 

New samples were combined with sequences available 
in GenBank, including our sequences previously submit- 
ted (Balakirev and Rozhnov 2010; Balakirev et al. 2011, 
2012, 2014; Rowe et al. 2008; Pages et al. 2010; Zhang 
et al. 2016). In total, 79 specimens were investigated (47 
C. langbianis, 15 C. thomasi and 17 C. chiropus). The 
geographic scope of the survey included 23 localities 
(Fig. 1, Table 1; see also Suppl. material 1: Table S1) 
scattered over China, Laos, Vietnam and Cambodia and 
constitute the known species distribution. Most of these 
sample specimens were collected personally by the au- 
thors (AVA, AEB). 


DNA extraction 


Small pieces of liver or muscle tissue were sampled in the 
field and stored in 96% ethanol. Total genomic DNA was 
extracted using a routine phenol/chloroform/proteinase K 
protocol (Kocher et al. 1989; Sambrook et al. 1989). The 
DNA was further purified either by a DNA Purification 
Kit (Fermentas, Thermo Fisher Scientific Inc., Pittsburgh, 
PA, USA) or by direct ethanol precipitation. We targeted 
five genes that were previously used for the phylogenetic 
analysis of Chiromyscus and were available for compar- 
ative analyses in GenBank. These genes included a com- 
plete Cytochrome b (Cyt b, 1140 bp); the 5’-proximal 
portion (680 bp) of subunit I of the Cytochrome Oxidase 
(COD), which 1s generally used for species diagnoses and 
for DNA barcoding for a number of mammals (Hebert et 
al. 2003); a portion of the first exon of the Interphotore- 
ceptor Retinoid Binding Protein gene (/RBP, also known 
as Rbp3, up to 1233 bp); the first exon of the Recombina- 
tion Activation Factor gene (RAG/, 1244 bp); and a por- 
tion of exon 10 of the Growth Hormone Receptor gene 
(GHAR, 815 bp). 


PCR amplification and sequencing 


Cyt b was amplified using H15915R, CytbRglu (Kocher 
et al. 1989; Irvin et al. 1991) and CytbRCb9H primers 
(Robins et al. 2007). The CO/ gene was amplified us- 
ing the universal conservative primers BatL 5310 and 
R6036R (Kocher et al. 1989; Irwin et al. 1991). The 
following universal PCR protocol was used to amplify 
mtDNA fragments: initial denaturation for 1 min 30 sec 
at 95 °C, 35 cycles of denaturation for 30 sec at 95 °C, 
annealing for 1 min at 52 °C and elongation for 30 sec 
at 72 °C, followed by terminal elongation for 2 min at 
72 °C. The PCR was performed in a 30-50 ul volume 
that contained 2.5—3 ul 10 x standard PCR buffer, 50 mM 
of each dNTP, 2 mM MgCl, 10 pmol of each primer, 
1 unit of Taq DNA polymerase (Fermentas, Thermo 
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Figure 1. Scattering of Chiromyscus spp. genetic samples in eastern Indochina. Type localities of C. chiropus, C. langbianis, and 


C. thomasi are shown by black circles. 


Fisher Scientific Inc., Pittsburgh, PA, USA) and 0.5 ul 
(25-50 ng) of total DNA template per tube. The reaction 
was performed using a Tercik (DNK-Tehnologia, Novo- 
sibirsk, Russia) thermocycler. JZRBP (1000-1610 bp in 
length) was amplified using the IRBP125f, IRBP1435r, 
IRBP1125r and IRBP1801r primers according to Stan- 
hope et al. (1992). A nested PCR technique was applied 
to amplify GHR in accordance with Jansa et al. (2009). 
An approximately 1.0-kb portion of exon 10 from the 
GHR gene was amplified using the primers GHRF1 and 
GHRendAlt. This PCR product was re-amplified using 
the nested GHRF1 primer paired with GHR750R and the 


GHRF50 primer paired with GHRendAlt. A 1244 bp por- 
tion of the RAG/ gene was obtained using primers S70 
and S71, as described by Steppan et al. (2004). 

The PCR products were purified using a DNA Purifica- 
tion Kit (Fermentas, Thermo Fisher Scientific Inc., Pitts- 
burgh, PA, USA). The resulting double-stranded DNA 
products were directly sequenced in both directions us- 
ing the Applied Biosystems 3130 Genetic Analyzer with 
the BigDye Terminator Cycle Sequencing Ready Reac- 
tion Kit (Applied Biosystems, Waltham, Massachusetts, 
USA). All obtained sequences were deposited in Gen- 
Bank (www.ncbi.nlm.nih.gov/Genbank; MK957014— 
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Table 1. List of geographical localities of Chiromyscus specimens used for genetic and morphological analyses. 


No. (Fig. 1) Species DNA lineage Locality Elevation (m asl) Latitude / Longitude 
1 C. langbianis N lineage China, Yunnan, Xishuangbanna 22.0°N, 100.8°E 
5 C. langbianis N lineage Vietnam, Tuyen Quang, Khong May 102 221soaN; 105, 339°F 
6 C. langbianis N lineage Vietnam, Vinh Phuc, Tam Dao 850 21.452°N, 105.636°E 
7 C. langbianis N lineage Vietnam, Lang Son, Huu Lien 230 21.661°N, 106.362°E 
8 C. langbianis N lineage Vietnam, Nghe An, Pu Hoat 840 19,.756°N, 104.796°E 
ile} C. langbianis N lineage Laos, Knhammouane TBSrN OS .30; 2 
13 C. langbianis N lineage Laos, Khammouane, Pha Deng 1.7.57 N; LObG23*E 
14 C. langbianis N lineage Vietnam, Quang Binh, Le Thuy, Sa Khia 156 17.068°N, 106.601°E 
ils: C. langbianis N lineage Vietnam, Kon Tum, Kon Plong 1030 14.722°N, 108.316°E 
16 C. langbianis N lineage Vietnam, Kon Tum, Kon Chu Rang 1020 14.505°N, 108.541°E 
17 C. langbianis N lineage Vietnam, Gia Lai, Kon Ka Kinh 900 14.203°N, 108.315°E 
19 C. langbianis S lineage Vietnam, Lam Dong, Bi Doup-Nui Ba 1400-1800 12.179°N, 108.679°E 
11 C. langbianis Hainan lineage China, Hainan, Jianfengling 18.74°N, 108.85°E 
12 C. langbianis Hainan lineage China, Hainan, Baoting 18.641°N, 109.775°E 
18 C. langbianis Cambodian lineage Cambodia, Kaoh Kong,Thmar Bang, Tatai Leu 11.961°N, 103.303°E 
3 C. thomasi N lineage Vietnam, Lao Cai, Bat Xat, Y Ty 1830 22.624°N, 103.629°E 
4 C. thomasi N lineage Vietnam, Son La, Muong Coi 547 21.343°N, 104.749°E 
5 C. thomasi N lineage Vietnam, Tuyen Quang, Khong May 102 2233SS NOS 3397 
2 C. thomasi S lineage Laos, Houay Sai, Houay Khot Station 20.267°N, 100.4°E 
9 C. thomasi S lineage Vietnam, Nghe An, Xoong Con 141 19.252°N, 104.318°E 
10 C. thomasi S lineage Vietnam, Nghe An, Pu Mat 200 18.957°N, 104.686°E 
13 C. thomasi S lineage Laos, Knhammouane, Pha Deng 17.57°N, 105.23°E 
tS C. thomasi S lineage Vietnam, Kon Tum, Kon Plong 1030 14.722°N, 108.316°E 
20 C. chiropus Vietnam, Lam Dong, Bao Loc 650 11.837°N, 107.64°E 
21 C. chiropus Vietnam, Dong Nai, Ma Da 75 11.381°N, 107.062°E 
22 C. chiropus Vietnam, Tay Ninh, Lo Go Xa Mat 1 5SS Ni 103,993°E 
23 C. chiropus Vietnam, Ba Ria-Vung Tau, Binh Chau 68 10.55°N, 107.483°E 


MK957137). Niviventer spp., Rattus norvegicus and Mus 
musculus were used as outgroups. 


Molecular data analyses 


Individual sequences were edited manually using BioEd- 
it v. 7.1.11 (Hall 1999) and aligned by Clustal W soft- 
ware incorporated into BioEdit and MEGA 6. The basic 
sequence parameter calculations and the best-fitting evo- 
lution models and inter- and intrapopulation divergence 
evaluations were performed using MEGA 6 (Tamura et al. 
2013). No pseudogenes were detected for the mitochon- 
drial genes. The optimal substitution models and their pa- 
rameters are summarised in Suppl. material 1: Table S2. 
Genetic distances (d) between groups under Tamuta-Nei 
gamma distributed invariant sites including (TN93+G+I) 
or Tamura 3-parameter (T3P) models (Tamura et al. 
2004), (depending on the best model determined) were 
calculated, based on the Cytb and COI genes in MEGA 6. 
Bayesian phylogenetic trees were inferred using MrBayes 
v.3.2. (Huelsenbeck and Ronquist 2001; Ronquist and 
Huelsenbeck 2003), two MCMCs for four chains with the 
default heating value and with a burn-in parameter equal 
to 25% of the initial number of runs. We applied 6x10° 
generations for the Cyt b dataset, 210° for the COJ, 4x10° 
for the RAG/ datasets and 5x10° for both GHR and IRBP 
datasets until the average standard deviation of split fre- 
quencies dropped below the level of 0.0025 after the runs 
for all datasets investigated. We used a flat Dirichlet prior 
for the relative nucleotide frequencies and for the relative 
rate parameters, a discrete uniform prior for the topolo- 
gies and an exponential distribution for the gamma shape 
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parameter and all branch lengths. The gamma shape pa- 
rameters for Bayesian Inference were evaluated directly 
from a general dataset by MrBayes v.3.2. A burn-in peri- 
od of one million generations was determined graphically 
using TRACER v.1.4 (Rambaut and Drummond 2007) to 
ensure convergence. Consensus trees were built from the 
last 25% of trees obtained (15 x10*, 15 x10*, 50 x10°, 10 
x10° and 12.5 x10? trees for Cyt b, COL, RAGI, GHR and 
IRBP, respectively) during the MCMC procedure by Mr- 
Bayes v.3.2. The five individual genes were concatenated 
using the software SequenceMatrix v1.7.6 (Vaidya et al. 
2011) to create a master alignment of 5,199 bp total (5,208 
bp including three triplet insertions in Mus musculus GHR 
gene). A restricted dataset, including all species includ- 
ed in this study and consisting of samples with a com- 
plete data matrix, were used for concatenated sequences 
analyses. A total of 5x10° generations was applied during 
the MCMC procedure by MrBayes v.3.2. for concatenat- 
ed alignment until the average standard deviation value 
dropped to 0.0077. TREEROT v.3 (Sorenson and Franzo- 
sa 2007) was used to examine tree-bisection-reconnection 
branch-swapping (PBS) to assess the contribution of each 
data partition in the combined analysis (Baker and De- 
Salle 1997). This analysis was performed to test the sus- 
tainability of the primary internal nodes for the different 
gene analyses. The robustness of the trees was assessed 
by posterior probabilities (PP). Trees were visualised and 
prepared by FigTree v.1.4.3 (Rambaut 2012). 

Divergence time approximation was performed by 
Mega X (Kumar et al. 2018), a time tree inferred using 
the Reltime method (Tamura et al. 2012; 2018) and the 
General Time Reversible model and branch lengths eval- 
uated by MrBayes v.3.2. for the concatenated dataset. The 
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Figure 2. The phylogenetic tree (Cyt b, Bayesian inference) for Chiromyscus genetic lineages radiation. The posterior probability 
values are presented at the nodes, and the branch lengths (scale bar at the bottom) are indicated above the nodes. The sample labels 
and locality numeration are indicated as in Fig. 1 and Suppl. material 1: Table S1. 


timetree was computed using one calibration constraint 
chosen as the divergence event between Mus and Rattus 
genera which are known to happen within 12.3—11.0 Mya 
(95% CI) (Benton and Donoghue 2007 with correction 
of Kimura et al. 2015). Discrete Gamma distribution was 
used to model evolutionary rate differences amongst sites 
(five categories (+G, parameter = 0.9185)). The rate vari- 
ation model allowed for some sites to be evolutionarily 
invariable ([+/], 54.21% sites). The analysis involved all 
82 nucleotide sequences and 5208 nucleotide positions. 


Morphological data analyses 


In total, 63 intact skulls of adult Chiromyscus (15 C. chiro- 
pus, 35 C. langbianis and 13 C. thomasi) obtained from 19 
genetically-investigated localities in Vietnam (see Table 1 
and Suppl. material 1: Tables S1 for references) were mea- 
sured for morphometric comparison and analyses. Age was 
assessed by tooth wear and closure of cranial sutures. Due 
to the limited sampling, sexual differences were not espe- 
cially tested; the possible sexual bias was compensated for 
by equalisation of representatives of different sexes in the 
sample. The sex ratio did not exceed 15 percentage points. 

Twenty measurements were taken from each skull by 
means of digital calipers to the nearest 0.01 mm: greatest 
length of skull (ONL), braincase breadth (BBC), brain- 
case height (HBC), zygomatic breadth (ZB), interorbital 
breadth (IB), length of rostrum (LR), breadth of rostrum 


(BR), breadth of zygomatic plate (BZP), diastema length 
(LD), length of foramina incisive (LIF), breadth of foram- 
ina incisive (BIF), length of bony palate (LBP), breadth 
across the palatal bridge at the level of the first molar 
(BBP), distance from the anterior edge of the premaxillary 
to the posterior edge of the palatine (= postpalatal length, 
PPL), breadth of the mesopterygoid fossa (BMF), length 
of the bulla (LB), upper molar row length (CLM 1-3), first 
upper molar breadth (BM1), first lower molar breadth 
(Bm1) and lower molar row length (CLm1-3). The cranial 
measurements followed Musser et al. (2006) and Balaki- 
rev et al. (2011), Suppl. material 1: Fig. S1. The measure- 
ments dataset is available from AEB by request. 

Principal components analysis (PCA) and canonical 
discriminant function analysis (DFA) were used to eval- 
uate “distinctiveness” amongst the samples. A one-way 
analysis of variance (ANOVA) was performed to test the 
differences amongst groups on all cranial variables. The 
software programme Statistica 8.0 (StatSoft Inc., Tulsa, 
OK, USA) was used for all analytical procedures. 


Results 


Phylogenetic subdivision and relationships 

The most representative tree was constructed for 66 
Cyt 6b sequences. The trees were well supported (PP = 
1) (Fig. 2). Different geographical populations of 
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C. chiropus formed a single homogeneous cluster, while 
C. thomasi and C. langbianis populations were represent- 
ed by a few geographically-segregated clusters. C. thom- 
asi split into two clusters, which we conventionally called 
the Northern and Southern phylogroups. The haplotypes 
forming the Northern cluster were distributed locally 
over the north-western part of Vietnam, whereas South- 
ern haplotypes spread substantially more widely, stretch- 
ing to the south to the Tay Nguyen Plateau (also known 
as the Central Highlands) and westwards to the extreme 
northwest of Laos. Chiromyscus langbianis formed three 
coherent geographically-distributed clusters. The first 
from the Dalat Plateau, the terra typica of this species; 
named hereafter as the Southern lineage. Samples origi- 
nating from continental Indochina, southern Yunnan and 
north-eastern Vietnam in the north to the Tay Nguyen Pla- 
teau in the south formed a second cluster, the Northern 
lineage. Another cluster, which was sister to the Northern 
lineage, contained samples from Hainan Island (Figure 
2). The genetic divergence within and between these 
groups is shown in Tables 2, 3. The estimated genetic dis- 
tances (d) for intraspecific C. langbianis and C. thomasi 
phylogroups were not very high and fell within the limits 
of 0.036 and 0.063. 

Another tree constructed using the mitochondrial CO/ 
gene of 43 samples revealed another additional specific 
phylogenetic lineage (Suppl. material 1: Figure S2), with 
divergence levels as high as those for the groups described 
above. (Table 3). It was formed by a single sample from 
the Cardamom Mountains, Cambodia. Unfortunately, this 
was the only sample of this group and, so far, there are 
no data about the distribution of this genetic lineage in 
southern Indochina. The clustering pattern of C. thomasi 
samples was similar to that obtained using Cyt b, with 
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C. chiropus demonstrating two reliable subclusters, com- 
bining the animals from the Dalat Plateau foothills (Lam 
Dong Province) and lowland populations (Tay Ninh, 
Dong Nai and Ba Ria-Vung Tau Provinces), respectively. 
Thus, mitochondrial genes supported C. /angbianis split 
into four geographically-distributed phylogroups, where- 
as C. thomasi and C. chiropus formed two distinctive 
phylogroups each. 

Phylogenetic reconstructions, based on nuclear genes, 
did not allow us to clarify the relationships and the tax- 
onomic rank of the distinctive phylogroups identified. 
Thus, only species-level clusters were reliably traced 
by the RAG/ gene tree constructed for the 26 available 
samples (Suppl. material 1: Fig. S3). The overall level 
of divergence was low and genetic distances did not ex- 
ceed 0.01. The same species-level groups were identified 
by the GHR gene (Suppl. material 1: Fig. S4), of which 
we included 52 samples. Within C. /angbianis, complex 
soft polytomy without notable geographic segregation 
was traced, whereas within C. thomasi, two clusters cor- 
responding to the mitochondrial phylogenetic lineages 
mentioned above were clearly demonstrated. At the same 
time, the considerable length of the branches was apparent 
for C. thomasi and for some specimens of C. langbianis 
from the Dalat Plateau. These branches were significantly 
longer than those characteristic of C. chiropus and most 
of the C. langbianis samples, a trait that may indicate a 
special pattern of its evolutionary history and, in partic- 
ular, the longer evolutionary age of these populations. 
Species-level clusters may also be traced in the /RBP 
gene tree, of which we had 49 samples (Suppl. material 
1: Fig. $5). Within the C. /angbianis cluster, no geograph- 
ical segregation was traced, which may indicate incom- 
plete sorting of lineages; on the other hand, C. thomasi 


Table 2. Genetic distances (d, TN93+G+I, gamma = 1.48) for geographic populations and species of Chiromyscus as calculated 
based on the Cyt b gene sequence (1140 bp). Standard error (S.E.) estimates are shown above the diagonal. 


C.llangbianis  C.langbianis —C_langbianis C.chiropus C.thomasi C.thomasi 
(Hainan) (Northern) (Southern) (Northern) (Southern) 
between group distances within groups distances 
d (TN93+G+I, Tamura-Nei) S.E. 

C.langbianis (Hainan) 0.005 0.008 0.012 0.015 0.016 0.0061 0.0013 
C.langbianis (Northern) 0.036 0.007 0.011 0.015 0.015 0.0097 0.0015 
C.langbianis (Southern) 0.063 0.052 0.013 0.016 0.016 0.0080 0.0017 
C.chiropus 0.129 0.123 0.130 0.014 0.013 0.0074 0.0015 
C.thomasi (Northern) 0.179 0.188 0.187 0.166 0.008 0.0039 0.0013 
C.thomasi (Southern) 0.186 0.191 0.192 0.162 0.063 0.0071 0.0016 


Table 3. Genetic distances (d; T3P, T92+I) for geographic populations and species of Chiromyscus as calculated based on the CO/ 
gene sequence (680 bp). Standard error (S.E.) estimates are shown above the diagonal. 


C.langbianis C.langbianis C.langbianis C.chiropus  C.chiropus — C.thomasi C.thomasi 
(Northern) (Southern) (Cambodia) (Lam_Dong) (others) (Northern) (Southern) 
between group distances within groups distances 
d(T3P: GTR) S.E. 

C.langbianis (Northern) 0.0106 0.0100 0.0176 0.0243 0.0289 0.0305 0.0045 0.0020 
C.langbianis (Southern) 0.035 0.0128 0.0183 0.0243 0.0296 0.0278 0.0022 0.0021 
C.langbianis (Cambodia) 0.033 0.044 0.0190 0.0242 0.0264 0.0292 
C.chiropus (Lam_Dong) 0.083 0.089 0.093 0.0121 0.0262 0.0277 0.0050 0.0024 
C.chiropus (others) 0.124 0.130 0.127 0.040 0.0294 0.0307 0.0022 0.0021 
C.thomasi (Northern) 0.156 0.161 0.138 0.155 0.172 0.0131 0.0110 0.0036 
C.thomasi (Southern) 0.163 0.153 0.154 0.160 0.178 0.049 0.0037 0.0026 
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Figure 3. A. The phylogenetic time tree (Cyt b/COY/RAGI/GHR/IRBP genes, concatenated analyses Bayesian inference) for Chiro- 
myscus genetic lineages radiation. The posterior probability values and average divergence time (Mya, in brackets) are presented at 
the nodes. Branches lengths are indicated above the branches. B. The position of Chiromyscus among most closely relative groups 
of rodents of SE Asia, marked by arrow (see Pages et al. 2016 for details). Footnote: The sample labels and locality numeration are 


indicated as in Fig. 1 and Suppl. material 1: Table S1. 


showed deep trichotomy. The samples corresponding to 
the Southern cluster of mitochondrial lineages were repre- 
sented here by two independent branches, which formed 
populations from the Tay Nguyen Plateau and populations 
distributed further to the north. Similar to the GHR gene 
tree (Suppl. material 1: Fig. S4), the inequality of the phy- 
logenetic branch lengths should garner attention. Howev- 
er, In contrast to the GHR gene, all samples of C. /angbi- 
anis from the Dalat Plateau recovered longer branches. In 
general, it can be concluded that C. /angbianis and C. chi- 
ropus showed significant homogeneity in nuclear genes, 
whereas C. thomasi had the same pattern of nuclear gene 
variation as traced by mitochondrial markers. 


In addition to low support levels for many nuclear 
gene clades, tree-bisection-reconnection branch-swap- 
ping (PBS) analysis indicates an existence of conflict- 
ing phylogenetic signals, especially for segments within 
C. langbianis. In general, the low posterior probability 
values for internal branches and the conflicting phylo- 
genetic signals in many lineages can be explained by a 
significantly slower evolution rate of nuclear genes (gen- 
erally weak phylogenetic signal) and incomplete lineage 
sorting that may be the result of symplesiomorphy. The 
tree which constructed the concatenated sequence (Fig. 3) 
is consistent with nuclear gene trees, but posterior proba- 
bilities values for some internal nodes are lower, mainly 
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Figure 4. Results of the multivariate analyses of Chiromyscus spp. from eastern Indochina. A. Ungrouped morphometric separation 
(PCA analysis); the data were drawn from 20 craniodental measurements. B. Grouped morphometric separation (DFA analysis) 


drawn from the same specimens and measurements. 


inside C. /angbianis. The Hainan cluster is not monophy- 
letic, as one of the samples was recovered in the Northern 
continental cluster. 


Morphological analyses 


The descriptive statistics of the craniodental measure- 
ments for phylogenetic lineages of C. /angbianis (two 
of four phylogroups discovered were available) and 
C. thomasi that were identified by the abovementioned 
analyses are summarised in Suppl. material 1: Tables S3, 
S4. Craniodental measurements for C. chiropus are given 
in Suppl. material 1: Table S5. As revealed by the inter- 
group F-test, the populations under study demonstrated 
notable peculiarities of cranial morphology. The samples 
were significantly different (p < 0.05 or lower) from each 
other in 18 and 7 cranial characters for C. Jangbianis and 
C. thomasi, respectively. 

In a principal components analysis (PCA) drawing on 
20 craniodental measurements, the first two axes captured 
60.9% (mainly reflecting general size) and 6.4% of the 
total variation, respectively. ONL, ZB, IB, LD, PPL and 
CLM1-3 were the six measurements that had the high- 
est correlations with PC 1 (Suppl. material 1: Table S6). 
In the PCA of skull measurements, all three species of 
Chiromyscus overlapped and C. /angbianis showed the 
largest range of variation amongst these species (Fig. 4A). 
Discriminant function analysis (DFA), which drew on the 
same variables, provided another means of illuminating 
the morphometric distinctions and the first two axes cap- 
tured 53.6% and 25.4% of the variation (Suppl. material 
1: Table S3). The DFA yielded moderate to high discrimi- 
nation amongst all species and genetic lineages (Fig. 4B). 
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Discussion 


Taxonomic implications 


The concordance of morphological and genetic traits and 
a good separation of samples in 3D factor space indicate 
the morphological specificity of the studied populations. 
On the other hand, the concordant pattern of morphologi- 
cal, genetic and clear geographic subdivision of the mito- 
chondrial phylogroups allow us to question the taxonom- 
ic status of these populations; in particular, they allow us 
to attribute the observed genetic lineages to distinct taxa. 

The Northern genetic lineage of C. Jangbianis must be 
undoubtedly assigned to subspecies C. /. indosinicus Os- 
good, 1932. This taxon was described as Rattus indosini- 
cus by Osgood (1932) from Sapa in northern Vietnam, Lao 
Cai Province. The Northern genetic lineage of C. thom- 
asi has to be attributed to the nominotypical subspecies 
C. t. thomasi (this phylogroup includes the holotype of 
C. thomasi). The appropriate name for the Southern ge- 
netic lineage of C. thomasi is debatable. It may be associ- 
ated with Rattus indosinicus vientianensis Bourret, 1942, 
described from the surroundings of Vientiane, Laos. How- 
ever, Musser (1973) treated vientianensis as a younger 
synonym of /angbianis and, in our previous survey where 
the most recent genus revision has been made (Balakirev 
et al. 2014), we also supposed that nomen vientianensis 
should be associated with C. /angbianis. Unfortunately, 
we have no specimens from the neighbouring Vientiane 
and we cannot identify which of the species is distributed 
there. The holotype of vientianensis is unavailable. 

The genetic distances between the two phylogroups of 
C. thomasi correspond minimally to the subspecific level 
(Baker 2006). However, despite their considerable ages, 
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Table 4. Estimated time to most recent common ancestor (Mya) for Chiromyscus based on Reltime method and the 
General Time Reversible model. A discrete Gamma distribution was used to model evolutionary rate differences among 
sites (5 categories (+G, parameter = 0.9185); The rate variation model allowed for some sites to be evolutionarily in- 


variable ({+/], 54.21% sites). The time tree was computed using | calibration constraints. 


Calibration Clade A Clade B Clade C Clade D Clade E Clade F Clade G 
Mus/Rattus Chiromyscus/ C. thomasi. C.chiropus _Northern/Southern Cambodian lineage Southern lineage Hainan lineage 
divergence Niviventercommon divergence’ divergence lineages of C. thomasi of C. langbianis of C. langbianis _ of C. langbianis 
point ancestor point point divergence point divergence point divergence point divergence point 
Mean 11.65 4.85 4.09 2.70 1.19 1.248 0.950 0.621 
95% CI lower 11.0 2.89 1.78 0.73 0.545 0.041 0.032 0.020 
95% CI upper 12.3 7.02 5.90 3.58 3.49 2.88 2.01 1.39 


these groups have no visually remarkable morphological 
differences (see Suppl. material 1: Table S3), so we can- 
not give here formal description of a new subspecies and 
assert only that the southern populations of C. thomasi 
belong to a unique monophyletic genetic lineage. 


Phylogeography and recent fauna formation 


Tree rats are usually confined to forest environments and 
their dispersal is restricted by the forest edge. They usual- 
ly do not spread beyond these limits and never cross wide 
deforested areas as they do not feel confident on open 
ground surface (Musser 1981; Corbet and Hill 1992; our 
personal observations). The main kind of natural event 
that contributed to the dispersal, population segregation 
and speciation of mammalian fauna is, most probably, 
repeated disjunction-reconnection events of natural pop- 
ulations that were associated with areas covered by tropi- 
cal forest during the late Miocene-Holocene (Hall 1998). 
The distribution pattern and genetic diversification of 
the genetic lineages revealed in the genus Chiromyscus 
shed light on the natural history and range formation of 
these species. Based on genetic data, the population of 
C. langbianis inhabiting the Dalat Plateau is apparently 
older than other recent continental populations. It can be 
assumed that the Dalat Plateau served as its main refu- 
gium during the Holocene climate oscillations (see also 
Meschersky et al. 2016; Abramov et al. 2017 for another 
rodent species). Judging by observed genetic distances 
and the homogeneity of the Northern genetic lineage of 
C. langbianis, its expansion beyond the Plateau occurred 
fairly quickly. Moreover, there is a reason to suppose 
multiple refugia, as supported by the fact that the most 
northern Hainanese cluster occupies a basal position in 
relation to the continental lineages. Noteworthy, to date, 
the Hainanese haplogroup may occur amongst continental 
populations. This may indicate both incomplete lineage 
sorting in this pair of clusters and an ancient hybridisation 
event between insular and continental populations. Multi- 
ple reconnection events that occurred during the Pleisto- 
cene make the second scenario possible. This gives some 
reason to believe that colonisation of the Island initiat- 
ed from a different population, not the one that inhabits 
the Dalat highlands. Instead, it probably originated from 
an additional northern refugium. The colonisation of 


Hainan by C. /angbianis might have happened simulta- 
neously with those of other Muridae (Niviventer and Rat- 
tus), Which are currently represented by distinct insular 
populations (Pan et al. 2007; Li et al. 2008; Smith and 
Xie 2008). This event could be dated back to the Late 
Miocene. According to Voris et al. (2000), Hainan had 
been connected to the mainland when the sea level was 
120-75 m below the current level, which has happened 
many times, with the longest connections occurring at 
approximately 0.25, 0.15 and 0.017 Mya. However, judg- 
ing by the estimated time of species level genetic lineage 
divergence (over | Mya, Fig. 3, Table 4), all of them were 
formed much earlier than these dates and cannot be asso- 
ciated with recent insularisation. It should be noticed that 
estimates, evaluated for divergence time for Muridae, are 
slightly higher than proposed earlier (Rowe et al. 2011; 
Pages et al. 2016); however, the genus Chiromyscus was 
represented there by only a single individual. Our finding 
provides evidence in support of more complex patterns of 
its evolutionary history. In any case, even if our estimates 
are closer to the higher limit of generic age determined 
earlier (Fabre et al. 2013; Pages et al. 2016), these timings 
for group split points are significantly older than the last 
events of the Hainan-Mainland reconnection. The latter 
supports the hypothesis of their formation in the con- 
tinental refugia during the Late Miocene. On the other 
hand, the occurrence of another original genetic lineage 
in southern Cambodia, which is an even more ancient 
separation than the Hainanese, indicates that there may 
have been several insularisation and resettlement events 
and that “distribution waves” originated from the Dalat 
and any other refugia during the Pleistocene. 

The split of C. thomasi into the Northern and South- 
ern phylogroups happened before the split of the corre- 
sponding C. /angbianis phylogroups and apparently is 
associated with antecedent global natural factor fluctu- 
ations. However, the recent distribution pattern of these 
Species indicates that their natural history differs signifi- 
cantly amongst the populations that originate from differ- 
ent dispersion centres/refugia. As far as can be traced by 
the data available, C. thomasi does not reach the Dalat 
Plateau and more southern regions inhabited by C. chi- 
ropus and the Southern lineage of C. /angbianis. At the 
same time, C. thomasi (both Northern and Southern phy- 
logroups) appears to be distributed sympatrically with the 
Northern phylogroup of C. /angbianis in most of eastern 
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Indochina. This indicates that their possible migration 
routes alongside the Annamite Range occurred in two op- 
posite directions, with C. thomasi moving northwards and 
C. langbianis moving southwards. The fact that C. thom- 
asi did not participate in mammal fauna formation on 
Hainan Island supports the recent natural area expansion 
of C. langbianis and are probably explained by ecological 
factors. Namely, these phenomena may reflect the eco- 
logical preferences of this species. C. thomasi is known 
to be more strictly associated with mountain forest forma- 
tions than C. langbianis, showing greater habitat versatil- 
ity, which apparently allowed the latter to spread much 
further eastwards along the plains of eastern Indochina. 
On the other hand, the significant genetic homogeneity of 
C. chiropus, which inhabits forest formations everywhere 
in the extreme south of Indochina and its basal position in 
relation to the genetic lineages of C. Jangbianis, may in- 
dicate that these recent populations diverged significantly 
earlier. This finding also indicates that forest refugia re- 
mained at the southern part of Indochina throughout the 
Holocene and even earlier. They could be associated not 
only with the Dalat Plateau, but also with the Cardamom 
Mountains, Bolaven Plateau and probably also with some 
of the offshore islands on the shelf of the Gulf of Siam. 

The distribution pattern of Chiromyscus species in the 
region also raises the problem of the initial intrusion and 
distribution of C. chiropus in eastern Indochina. The terra 
typica for this species is the Karen Mountains in eastern 
Myanmar, where the species inhabit mountainous forests. 
As we pointed out earlier (Balakirev et al. 2014), Burmese 
specimens do not demonstrate noticeable morphological 
differences from the southern Vietnamese populations and 
these are treated as conspecific. Unfortunately, there are 
still no genetic data from Burmese populations that would 
allow direct comparison of their genetic identity. Howev- 
er, the wide distribution of this species to the east in east- 
ern Indochina through the Yunnan and Annamite Ranges 
is hampered by the wide distribution of another species, 
namely, C. thomasi, which populates these regions. No 
cases of sympatry are currently documented, which may 
suggest a competitive exclusion in this pair of species. At 
the same time, the existence of a direct connection be- 
tween the Malacca and southern Indochina in the Holo- 
cene by a forest corridor cannot be excluded. Based on 
data of Meijaard et al. (2003) on tropical forest persistence 
and the distribution of forest-dependent species on islands 
of the South China Sea and a forest connection between 
southern Indochina and Malacca, a southern expansion 
route is probable. Nevertheless, there are no records on 
the current distribution of this species in the lowland areas 
in central Indochina to the west from 105°E. 


Conclusions 
We show that the genetic distances between phylogroups 


of C. langbianis and C. thomasi correspond to the subspe- 
cific level at least. However, these phylogenetic groups do 
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not demonstrate obvious univocal diagnostic differences 
in cranial features suitable for species diagnoses without 
special statistical analysis. Our study shows that the recent 
phylogenetic structure of C. Jangbianis is the most recent 
within the genus and appears within several independent 
refugia that remained isolated throughout the Pleistocene. 
In turn, the phylogroups of C. thomasi are likely older than 
those of C. /angbianis. Environmental factors and species 
preferences followed recent natural ecological shifts which 
drove allopatry. However, C. chiropus demonstrates the 
greatest age; the ways of formation of the area of this spe- 
cies still remain obscure and are likely to be associated with 
changes in forest cover in Indochina and Malacca Penin- 
sula during the Pleistocene. The possibility of competitive 
interaction of these species in the process of formation of 
their recent natural areas also cannot be excluded. 
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